# TODO: Add comment
# 
# Author: broken
###############################################################################
library(raster)
library(base)
library(gstat)
library(RPostgreSQL)
library(hydroGOF)
library(stringr)
dr=dbDriver("PostgreSQL")

path= "C:/Users/forever/Desktop/NPP01/"

DB = "fimo"
HOST ="192.168.0.4"
PORT = 5432
USER = "rasdaman"
PASS = "rasdaman"

con=dbConnect(dr,dbname=DB,host=HOST,port=5432,user=USER,password=PASS)
timeQuery = "SELECT aqstime FROM apom.satresampviirs where filename like '%"
selectStr = "SELECT stationid,name,st_asText(location),avg(pm25) as avg_pm FROM apom.grdpmhiscem_data,public.grdstation WHERE grdpmhiscem_data.stationid = grdstation.id AND grdstation.id IN (1,2,4,5) AND  aqstime between '"
andStr="' and '"
groupStr="'  group by stationid,name,location having avg(pm25) > 0"

data_result= data.frame()
year_array=c(2012,2013,2014)
for (k in 1:3) {
	inputPath = paste(path,year_array[k],sep="")
	#inputPath = paste(path,2012,sep="")
	listFile = list.files(path = inputPath,pattern="_rg.tif$",full.names = FALSE,recursive=TRUE)
	listFile = str_replace(listFile,"_rg.tif","")
	listSize=length(listFile)
	
	listTime=c(1:listSize)
		
	uCor=uRMSE=uRE=uMFE=uMFB = c(1:listSize)
	totalRow=c(1:listSize)
	
	for(i in 1:listSize){
		fileName=listFile[i]
		
		ogFile = paste(inputPath,"/",fileName,"_rg.tif",sep="")
		ukFile = paste(inputPath,"/",fileName,"_uk.tif",sep="")
	
		ogRaster=raster(ogFile)
		ukRaster=raster(ukFile)
	
		search_query = paste(timeQuery,fileName,"%'",sep="")
		resultSet = dbSendQuery(con,search_query)
		record = fetch(resultSet,n=1)
		asqTime = substr(record[1,1],1,30)
		myTime = strptime(asqTime,format="%Y-%m-%d %H:%M:%S")
	
		#add 7h
		myTime = myTime+25200
		
		startTime = myTime-1800
		endTime = myTime+1800
		
		aqsMonth = format.Date(myTime,"%m")
		aqsYear = format.Date(myTime,"%Y")
		
		selectQuery=paste(selectStr,startTime,andStr,endTime,groupStr,sep="")
		pmRs=dbSendQuery(con,selectQuery)
		pmRecord=fetch(pmRs,n=-1)
		numberRow=nrow(pmRecord)
		if(numberRow>0){
			ogValueArray=c(1:numberRow)
			ukValueArray=c(1:numberRow)
			fileArray = c(1:numberRow)
			yearArray = c(1:numberRow)
			monthArray = c(1:numberRow)
			aqsArray = c(1:numberRow)
		
			for(j in 1:numberRow){
				locationStr=pmRecord[j,3]
				index=regexpr(" ",locationStr)[1]
				xLocation=substr(locationStr,7,index-1)
				yLocation=substr(locationStr,index+1,nchar(locationStr)-1)
				
				ogValue=extract(ogRaster,cbind(as.numeric(xLocation),as.numeric(yLocation)))
				ukValue=extract(ukRaster,cbind(as.numeric(xLocation),as.numeric(yLocation)),buffer=25000,fun=mean)
				
				ogValueArray[j]=ogValue
				ukValueArray[j]=ukValue
				
				fileArray[j]=fileName
				yearArray[j]=aqsYear
				monthArray[j]=aqsMonth
				aqsArray[j]=myTime
						
						
			
			}
			pmRecord$ogColumn=ogValueArray
			pmRecord$ukColumn=ukValueArray
			
			pmRecord$fileName=fileArray
			pmRecord$year=yearArray
			pmRecord$month=monthArray
			pmRecord$aqstime=aqsArray
			
			
			data_result = rbind(data_result,data.frame(pmRecord))
			print(paste("Complete ",fileName,sep=" "))
		
		}else{
			print(paste(fileName,"have no data ",sep=" "))
		}
		
	}
	
	
}

totalRow=nrow(data_result)

u_cor=cor(data_result$avg_pm,data_result$ukColumn)*cor(data_result$avg_pm,data_result$ukColumn)
u_rmse=sqrt(sum((data_result$ukColumn-data_result$avg_pm)^2)/totalRow)
u_re=(sum(abs(data_result$ukColumn-data_result$avg_pm)/data_result$avg_pm)/totalRow)*100

u_mfb = (2*sum((data_result$ukColumn-data_result$avg_pm)/(data_result$avg_pm+data_result$ukColumn)*100))/totalRow
u_mfe = (2*sum(abs(data_result$ukColumn-data_result$avg_pm)/(data_result$ukColumn+data_result$avg_pm)*100))/totalRow


outFile=paste(path,"npp_avg_year_4np.csv",sep="")
write.csv(data_result,file=outFile)

u_cor
u_rmse
u_re

print("Finish All")




